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Abstract. We study the density of states in monolayer and bilayer graphene in the 
presence of a random potential that breaks sublattice symmetries. While a uniform 
symmetry-breaking potential opens a uniform gap, a random symmetry-breaking 
potential also creates tails in the density of states. The latter can close the gap again, 
preventing the system to become an insulator. However, for a sufficiently large gap 
the tails contain localized states with nonzero density of states. These localized states 
allow the system to conduct at nonzero temperature via variable-range hopping. This 
result is in agreement with recent experimental observations in graphane by Elias et 
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1. Introduction 

Graphene is a single sheet of carbon atoms that is forming a honeycomb lattice. A 
graphene monolayer as well as a stack of two graphene sheets (i.e. a graphene bilayer) 
are semimetals with remarkably good conducting properties [H EJ [3]. These materials 
have been experimentally realized with external gates, which allow for a continuous 
change in the charge-carrier density. There exists a non-zero minimal conductivity at 
the charge neutrality point. Its value is very robust and almost unaffected by disorder 
or thermal fluctuations [3], HI [6] . 

Many potential applications of graphene require an electronic gap to switch between 
conducting and insulating states. A successful step in this direction has been achieved by 
recent experiments with hydrogenated graphene (graphane) [7J and with gated bilayer 
graphene [El El HQ] • These experiments take advantage of the fact that the breaking of 
a discrete symmetry of the lattice system opens a gap in the electronic spectrum at the 
Fermi energy. In the case of monolayer graphene (MLG), a staggered potential that 
depends on the sublattice of the honeycomb lattice plays the role of such symmetry- 
breaking potential (SBP). For bilayer graphene (BLG) a gate potential that distinguishes 
between the two graphene layers plays a similar role. 

With these opportunities one enters a new field in graphene, where one can switch 
between conducting and insulating regimes of a two-dimensional material, either by a 
chemical process (e.g. oxidation or hydrogenation) or by applying an external electric 
field HI]. 

The opening of a gap can be observed experimentally either by a direct measurement 
of the density of states (e.g., by scanning tunneling microscopy [12]) or indirectly by 
measuring transport properties. In the gapless case we observe a metallic conductivity 
a oc pD, where D is the diffusion coefficient (which is proportional to the scattering 
time) and p is the density of states (DOS). This gives typically a conductivity of the order 
of e 2 /h. The gapped case, on the other hand, has a strongly temperature-dependent 
conductivity due to thermal activation of charge carriers [13] 

a(T) = a e- T °' T (1) 

with some characteristic temperature scale T which depends on the underlying model. 
A different behavior was found experimentally in the insulating phase of graphane [TJ: 

a(T) « a e~^ 1/3 , (2) 

which is known as 2D variable-range hopping [H]. This behavior indicates the 
existence of well-separated localized states, even at the charge- neutrality point, where 
the parameter T depends on the DOS at the Fermi energy E F as T oc l/p(E F ). 

The experimental observation of a metal-insulator transition in graphane raises two 
questions: (i) what are the details that describe the opening of a gap and (ii) what is 
the DOS in the insulating phase? In this paper we will focus on the mechanism of the 
gap opening due to a SBP in MLG and BLG. It is crucial for our study that the SBP 
is not uniform in the realistic two-dimensional material. One reason for the latter is 
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the fact that graphene is not flat but forms ripples [THJ [TjJJ [T7] . Another reason is the 
incomplete coverage of a graphene layer with hydrogen atoms in the case of graphane 
[7j. The spatially fluctuating SBP leads to interesting effects, including a second-order 
phase transition due to spontaneous breaking of a discrete symmetry and the formation 
of Lifshitz tails. 

2. Model 

Quasiparticles in MLG or in BLG are described in tight-binding approximation by a 
nearest-neighbor hopping Hamiltonian 

H = — t T y(? r c T i + ^ V r c\c r + h.c, (3) 

<r,r'> r 

where cj (c r ) are fermionic creation (annihilation) operators at lattice site r. The 
underlying lattice structure is either a honeycomb lattice (MLG) or two honeycomb 
lattices with Bernal stacking (BLG) [TTJ [18] . We have an intralayer hopping rate t and 
an interlayer hopping rate t± for BLG. There are different forms of the potential V r , 
depending on whether we consider MLG or BLG. Here we begin with potentials that are 
uniform on each sublattice, whereas random fluctuations are considered in subsection 

El 

2.1. MLG 

V r is a staggered potential with V r = m on sublattice A and V r = — m on sublattice 
B. This potential obviously breaks the sublattice symmetry of MLG. Such a staggered 
potential can be the result of chemical absorption of non-carbon atoms in MLG (e.g. 
oxygen or hydrogen [7]). A consequence of the symmetry breaking is the formation of 
a gap A g = m: The spectrum of MLG consists of two bands with dispersion 

E k = ±^m 2 + el (4) 

where 

e 2 k = t 2 [3 + 2 cos fci + 4 cos(Aa/2) cos( V3A&/2)] (5) 
for lattice spacing a — 1. 

2.2. BLG 

V r is a biased gate potential that is V r = m (V r = —m) on the upper (lower) graphene 
sheet. The potential in BLG has been realized as an external gate voltage, applied to 
the two layers of BLG [8]. The spectrum of BLG consists of four bands [11] with two 
low-energy bands 

EkH =±Je 2 + tl/2 + - ^1/4 + (t\ + 4m^ , (6) 
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where is the monolayer dispersion of Eq. (J5]), and two high-energy bands 

E+(m) = ±yj 4 + t\/2 + mH y^/4 + (t\ + 4m 2 )e 2 • (7) 

The spectrum of the low-energy bands has nodes for m = where E^(0) vanishes 
in a (k — K) 2 manner, where K is the position of the nodes, which are the same as 
those of a single layer. For small m <C t±, a mexican hat structure develops around 
k = K, with local extremum in the low-energy band at E^(m) = ±m, and a global 
minimum/maximum in the upper/lower low energy band at E^(m) = mt±/ \ft\ + 4m 2 . 

For small gating potential V r = ±m we can expand E^(m) under the square root 
near the nodes and get 

E-{m) ~ ±^[1 - 4e?A(ti + 4e 2 fc )-V>2 + S -( )a . (8) 

t± apparently reduces the gap. Very close to the nodes we can approximate the 
factor in front of m 2 by 1 and obtain an expression similar to the dispersion of MLG: 
±^m 2 + ^7(0) 2 . Here we notice the absence of the mexican hat structure 
in this approximation. The resulting spectra for MLG and BLG are shown in Fig. [TJ 




fcl - K 



Figure 1. The energy spectra of MLG (blue) and BLG (red) are shown, with and 
without a gap (dashed and solid line, respectively) for positive energies. Note the 
characteristic mexican hat structure of gapped BLG. 



2.3. Low-energy approximation 

The two bands in MLG and the two low-energy bands in BLG represent a spinor-1/2 
wave function. This allows us to expand the corresponding Hamiltonian in terms of 
Pauli matrices crj as 

H = hiai + h 2 o- 2 + ma 3 . (9) 
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Near each node the coefficients hj read in low-energy approximation [TJ] 

hj = iVj {MLG), h x = Vj - Vj, h 2 = 2ViV 2 (BLG) , (10) 
where (Vi, V2) is the 2D gradient. 

2.4- Random fluctuations 

In a realistic situation the potential V r is not uniform, neither in MLG nor in BLG, 
as discussed in the Introduction. As a result, electrons experience a randomly varying 
potential V r along each graphene sheet, and m in the Hamiltonian of Eq. (Q becomes a 
random variable in space as well. For BLG it is assumed that the gate voltage is adjusted 
at the charge-neutrality point such that in average m r is exactly antisymmetric with 
respect to the two layers: (mi) m = -(m2) m . 

At first glance, the Hamiltonian in Eq. (J3j) is a standard hopping Hamiltonian 
with random potential V r . This is a model frequently used to study the generic case of 
Anderson localization [20]. The dispersion, however, is special in the case of graphene 
due to the honeycomb lattice: at low energies it consists of two nodes (or valleys) K 
and K' [T7l [19] . It is assumed here that randomness scatters only at small momentum 
such that intervalley scattering, which requires large momentum at least near the nodal 
points (NP), is not relevant and can be treated as a perturbation. Then each valley 
contributes separately to the DOS, and the contribution of the two valleys to the DOS 
p is additive: p = px + Pk 1 - This allows us to consider the low-energy Hamiltonian in 
Eqs. ([9]), (fTOl) . even in the presence of randomness for each valley separately. Within 
this approximation the term m r is a random variable with mean value (m r ) m = m and 
variance ((m r — fh)(m r i — fh)) m = g5 r y. The following analytic calculations will be 
based entirely on the Hamiltonian of Eqs. (}9l.( lT0l) and the numerical calculations on 
the lattice Hamiltonian of Eq. ([3]). In particular, the average Hamiltonian (H) m can be 
diagonalized by Fourier transformation and is 

(H) m = k x a x + k 2 a 2 + rhcr 3 (11) 
for MLG with eigenvalues Ef. = ±\Jm 2 + k 2 . For BGL the average Hamiltonian is 

(H) m = {k\ - k\)a x + 2k ± k 2 a 2 + ma 3 (12) 
with eigenvalues = ±y/rh 2 + fc 4 . 

2.5. Symmetries 

Low-energy properties are controlled by the symmetry of the Hamiltonian and of the 
corresponding one-particle Green's function G(ie) = (H + ie)~ l . In the absence of 
sublattice-symmetry breaking (i.e. for m — 0), the Hamiltonian H = h\0\ + h 2 a 2 has a 
continuous chiral symmetry 



H -> e aa3 He aa3 = H 



(13) 
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with a continuous parameter a, since H anticommutes with er 3 . The term ma^ breaks 
the continuous chiral symmetry. However, the behavior under transposition hj = —hj 
for MLG and hj = hj for BLG in Eq. ( jTUl) provides a discrete symmetry: 

H -> -a n H T a n = H , (14) 

where n = 1 for MLG and n = 2 for BLG. This symmetry is broken for the one- 
particle Green's function G(ie) by the ie term. To see whether or not the symmetry is 
restored in the limit e — > 0, the difference of G(ze) and the transformed Green's function 
—cr n G T (ie)a n must be evaluated: 

G(ie) + a n G T (ze)a n = G(ie) - G(-ie) . (15) 

For the diagonal elements this is the DOS at the NP p(E = 0) = po in the limit e — > 0. 
Thus the order parameter for spontaneous symmetry breaking is po- According to the 
theory of phase transitions, the transition from a nonzero po (spontaneously broken 
symmetry) to po = (symmetric phase) is a second-order phase transition, and should 
be accompanied by a divergent correlation length at the transition point. Since our 
symmetry is discrete, such a phase transition can exists in d = 2 and should be of Ising 
type. A calculation, using the SCBA of po, gives indeed a second-order transition at the 
point where po vanishes with a divergent correlation length £ for the DOS fluctuations 

for m 2 ~ m 2 c with a finite coefficient £o [21]. Whether or not this transition is an 
artefact of the SCBA or represents a physical effect due to the appearence of two types 
of spectra (localized for vanishing SCBA-DOS and delocalized for nonzero SCBA-DOS) 
is not obvious here and requires further studies. 



2. 6. Density of states 

Our focus in the subsequent calculation is on the DOS of MLG and BLG. In the absence 
of disorder, the DOS of 2D Dirac fermions opens a gap A oc fh as soon as a nonzero 
term rh appears in the Hamiltonian of Eq. ([9]), since the low-energy dispersion is 
E k = ±Vrh 2 + k 2 for MLG and E k = ±^m 2 + k 4 for BLG, respectively (cf Fig. H}. 
Here we evaluate the DOS of MLG and BLG in the presence of a uniform gap. Given 
the energy spectrum, the DOS is defined as 

p{E) = Y,KE-E k ). (16) 

k 

By using the MLG dispersion, this reduces to 

p(E) = \E\Q(\E\-m), (17) 

where Q(x) is the Heaviside function. For BLG, this gives 

\E\ . . . 
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Figure 2. Density of states for a uniform symmetry-breaking potential for monolayer 
graphene and bilayer graphene is shown in the left panel. The density of states for a 
uniform symmetry-breaking potential for BLG is shown for several values of fj_. For 
small t±, the mexican hat structure influences the DOS by shifting the gap to lower 
values, and by developing a kink at E = m. 



which are shown in Fig. [21 By retaining the full low-energy spectrum for BLG, E k , the 
DOS can still be evaluated in closed form, with the result 

, {tl+Am2) , for m > \E\ > ,f 1 



p(E) = \E\x { / ; x , , (19) 

, , ( ±+ } , +1 for E > m. 

In the limit of t± ^> (E,m), this reduces to Eq. ffl8|) after dividing it by tj_, which was 
set to 1 in the low-energy approximation, and the DOS saturates to a constant value 
after the initial divergence. For finite t±, however, the Dirac nature of the spectrum 
appears again, and the high energy DOS increases linearly even for the BLG, similarly 
to the MLG case. For m — 0, and E <ti, this lengthy expression gives 

p(E « f x ) = t f. (20) 

An interesting question, from the theoretical as well as from the experimental point 
of view, appears here: What is the effect of random fluctuations around ml Previous 
calculations, based on the self-consistent Born approximation (SCBA), have revealed 
that those fluctuations can close the gap again, even for an average SBP term m ^ 
[22] . Only if m exceeds a critical value m c (which depends on the strength of the 
fluctuations), an open gap was found in these calculations. This describes a special 
transition from metallic to insulating behavior. In particular, the DOS at the Dirac 
point po vanishes with fa like a power law 

p (m) ~ ^Jm-m 2 c . (21) 
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The exponent 1/2 of the power law is probably an artefact of the SCBA, similar to the 
critical exponent in mean-field approximations. 



3. Self-consistent Born approximation 

The average one-particle Green's function can be calculated from the average 
Hamiltonian (H) m by employing the self-consistent Born approximation (SCBA) 

mm 



(G(ie)) m « ((H) m + ie- = G (iv, m.) . (22) 

The SCBA is also known as the self-consistent non-crossing approximation in the Kondo 
and superconducting community. The self-energy S is a 2 x 2 tensor due to the spinor 
structure of the quasiparticles: X = —(ir/ao + m s a 3 )/2. Scattering by the random SBP 
produces an imaginary part of the self-energy rj (i.e. a one-particle scattering rate) and 
a shift m s of the average SBP fh (i.e., fh — > m' = fh + m s ). £ is determined by the 
self-consistent equation 

£ = -ga 3 ((H) m + xe - 2S)"V 3 • (23) 

The symmetry in Eq. (1141) implies that with £ also 

(T n E(T n = -{irjao - m s a 3 )/2 (24) 

is a solution (i.e. m s — > — m s creates a second solution). 

The average DOS at the NP is proportional to the scattering rate: po = 
This reflects that scattering by the random SBP creates a nonzero DOS at the NP if 
r] > 0. 

Now we assume that the parameters rj and m s are uniform in space. Then Eq. (1231) 
can be written in terms of two equations, one for the one-particle scattering rate rj and 
another for the shift of the SBP m s , as 

7] = glr), m s = -rhgl/ (1 + gl) . (25) 

/ is a function of fh and 7] and also depends on the Hamiltonian. For MLG it reads with 
momentum cutoff A 

A 2 



I MLG — 7T~ l n 

Z7T 



1 + 



(26) 



and for BLG 

1 



l BLG 



rf + (m + m s ) 

(A ~ oo) . (27) 



4>/?7 2 + (m + m s ) 2 

A nonzero solution ri requires gl = 1 in the first part of Eq. (|25|) , such that m s = —fh / 2 
from the second part. Since the integrals I are monotonically decreasing functions for 
large fh, a real solution with gl = 1 exists only for \rh\ < m c . For both, MLG and BLG, 
the solutions read 

rj 1 = (m 2 c - m 2 )6(m 2 - rh 2 ) /4 , (28) 
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Figure 3. Schematic shape of the density of states: full curves are the bulk 
density of states for uniform symmetry-breaking potential, dotted curves represent the 
broadening by disorder. The broadened density of states can overlap inside the gap 
for to < m c (a) or not for to > m c (b), depending on the average symmetry-breaking 
potential to. m c is given in Eq. (129)) . 



where the model dependence enters only through the critical average SBP m c : 



m c is much bigger for BGL, a result which indicates that the effect of disorder is much 
stronger in BLG. This is also reflected by the scattering rate at rh = which is rj = m c /2. 

A central assumption of the SCBA is a uniform self-energy E. The imaginary 
part of E is the scattering rate t], created by the random fluctuations. Therefore, a 
uniform rj means that effectively random fluctuations are densely filling the lattice. If 
the distribution of the fluctuations is too dilute, however, there is no uniform nonzero 
solution of Eq. (|23|) . Nevertheless, a dilute distribution can still create a nonzero DOS, 
as we will discuss in the following: we study contributions to the DOS due to rare events, 
leading to Lifshitz tails. 

4. Lifshitz tails 

In the system with uniform SBP the gap can be destroyed locally by a local change of the 
SBP m ^m + 5m r due to the creation of a bound state. We start with a translational- 
invariant system and add Sm r on site r. To evaluate the corresponding DOS from the 
Green's function G = (H + ie + dmas)" 1 , using the Green's function Go = (H + ie)" 1 
with uniform m, we employ the lattice version of the Lippmann-Schwinger equation [26] 




MLG 
BLG 



(29) 




(30) 



(31) 
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In the case of MLG we have 

Go = [(E + , e y - H I £ {e _ iE)2 k +m2 + k2 dk (32) 

~ (Eao - ma 3 ) — log[l + A 2 /(m 2 - £ 2 )] + o(ie) = (g Q + zes)<7 + 03^33) 

47T 

(remark: the DOS of BLG has the same form.) Then the imaginary part of the Green's 
function reads 

/m[G(,)] = -( + * + ) (34) 

y o c) es (5f - 03 - H) y 

with 

= • (35) 

Thus the DOS is the sum of two Dirac delta peaks 

p r oc 5 ts (g + 03 + 5m r ) + 5 es (0 O - 03 - 5m r ) . (36) 

The Dirac delta peak appears with probability oc exp(— (00 ± 9^) 2 / 9) f° r a Gaussian 
distribution. This calculation can easily be generalized to 5m r on a set of several 
sites r (26]. Then the appearance of the several such Dirac delta peaks decreases 
exponentially. Moreover, these contributions are local and form localized states. For 
stronger fluctuations 5m r (i.e., for increasing 0) the localized states can start to overlap. 
This is a quantum analogue of classical percolation. 

The localized states in the Lifshitz tails can be taken into account by a 
generalization of the SCBA to non-uniform self-energies. The main idea is to search 
for space-dependent solutions S r of Eq. (1231) . In general, this is a diffult problem. 
However, we have found that this problem simplifies essentially when we study it in 
terms of a 1/fh expansion. Using a Gaussian distribution, this method gives Lifshitz 
tails of the form 



»<*> ^ bS^^ 1 ■ (37) 

5. Numerical approach 

To understand to behavior of random gap fluctuations in graphene, and also the 
limitations of the SCBA, we carried out extensive numerical simulations on the 
honeycomb lattice, allowing for various random gap fluctuations on top of a uniform 
gap m. These fluctuations are simulated by box and Gaussian distributions. From the 
SCBA, the emergence of a second-order phase transition at a critical mean m c is obvious 
for a given variance. This is best manifested in the behavior of the DOS, which stays 
finite for (m) < m c , and vanishes afterwards, and serves as an order parameter. Does 
this picture indeed survive, when higher order corrections in the fluctuations are taken 
into account? 
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To start with, we take a fix random mass configuration with a given variance and 
the honeycomb lattice (HCL) with the conventional hoppings (t), represented by H . 
Then, we take a separate Hamiltonian, responsible for the uniform, non-fluctuating gaps, 
denoted by H gap , and study the evolution of the eigenvalues of H + mH gap by varying 
m for a 600x600 lattice. By using Lanczos diagonalization, we focus our attention 
only to the 200 eigenvalues closest to the NP. Their evolution is shown in Fig. HI 
This supports the existence of a finite m c , but since it originates from a single random 
disorder configuration, rare events can alter the result. As a possible definition of the 
rigid gap, we also show the maximum of the energy level spacing for these eigenvalues 
as a function of m. As seen, it starts to increase abruptly at a given value of m, which 
can define m r . 




0.5 1 1.5 2 

m/t 



Figure 4. (Color online) The evolution of the 200 lowest eigenvalues is shown for a 
given random mass configuration with Gaussian distribution (with variance g) on a 
600x600 HCL, by varying the uniform gap. The red line denotes the maximum of the 
level spacing of these eigenvalues, a possible definition of the average gap. 

To investigate whether a finite critical m c survives, we take smaller systems and 
evaluate the averaged DOS directly from many disorder realizations. To achieve this, we 
take a 200x200 HCL, and evaluate the 200 closest eigenvalues to the NP, and count their 
number in a given small interval, AE (smaller than the maximal eigenvalues) around 
zero. This method was found to be efficient in studying other types of randomness [28J. 
We mention that large values of AE take contribution from higher energy states into 
account, while too small values are sensitive to the discrete lattice and consequently the 
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Figure 5. (Color online) The density of states at the NP is plotted for Gaussian 
random mass for a 200x200 HCL for g = 0.9 2 , 1, l.l 2 , 1.2 2 and 1.3 2 from bottom to 
top after 400 averages. The symbols denote the numerical data, solid lines are fits 
using a cxp(— bm c ). The inset shown the obtained exponents, c, as a function of g, 
which is close to 1.5. 



discrete eigenvalue structure of the Hamiltonian. For lattices containing a few 10 4 — 10 5 
sites, AE/t ~ 10~ 2 — 10 -4 are convenient. 

The resulting DOS is plotted in Figs. [5] and [6] for Gaussian (with variance g) and 
box distribution (within [— W..W], variance g = W 2 /3). This does not indicate a sharp 
threshold, but rather the development of long Lifshitz tails due to randomness, as we 
already predicted in the previous section. To analyze them, we fitted the numerical data 
by assuming exponential tails of the form 

p(0) = texp(-a - bm c ) (38) 

for a Gaussian and 

p(0) =teKp(-a-b/\m- W\ c ) (39) 

for a box distribution, as suggested by Ref . [29] . The obtained c values are visualized in 
the insets of Figs. [5]and[6j Given the good agreement, we believe that the DOS at the 
NP is made of states that are localized in a Lifshitz tail. We mention that these results 
are not sensitive to finite size scaling at these values of the disorder and uniform gap, 
only smaller systems (like the 30x30 HCL) require more averages (~ 10 4 ), whereas for 
larger ones (such as the 200x200 with 400 averages) fewer averages are sufficient. 

In Fig. [7J the energy dependent DOS is shown for Gaussian random mass with 
g = 1 and for several uniform gap values. With increasing m, the DOS dimishes rapidly 
at low energies, and develops a pseudogap. The logarithmic singularity at E — t is 
washed out for g — 1. We also show the inverse of the DOS, proportional to T , the 
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characteristic temperature scale of variable range hopping as a function of the carrier 
density (which is proportional to E 2 ). 



0.14 
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Figure 6. (Color online) The density of states at the NP is plotted for box distributed 
{[-W..W]) randomness for a 200x200 HCL for W = 1.5, 1.7 and 2 (g = W 2 /3) from 
bottom to top after 400 averages. The symbols denote the numerical data, solid lines 
are fits using aexp(— b/\m — W^| c ). The inset shown the obtained exponents, c, as a 
function of g. 



6. Discussion 

MLG and BLG consist of two bands that touch each other at two nodal points (or 
valleys). Near the nodes the spectrum of MLG is linear (Dirac-like) and the spectrum 
of BLG is quadratic. The application of a uniform SBP opens a gap in the DOS for both 
cases. For a random SPB, however, the situation is less obvious. First of all, it is clear 
that randomness leads to a broadening of the bands. If we have two separate bands 
due to a small uniform SPB, randomness can close the gap again due to broadening (cf. 
Fig. (2^,). The broadening of the bands depends on the strength of the fluctuations of 
the random SBP. In the case of a Gaussian distribution there are energy tails for all 
energies. 

Now we focus on the NP, i.e. we consider E = and po- Then we have two 
parameters in order to change the gap structure: the average SBP (m) = fa and the 
variance g. fa allows us to broaden the gap and g has the effect of closing it due to 
broadening of the two subbands. Previous calculations have shown that at the critical 
value m c (g) of Eq. (129]) the metallic behavior breaks down for fa > ra c (g) [22]. On 
the other hand, Gaussian randomness creates tails at all energies. Consequently, there 
are localized states for \fa\ > m c (g) at the NP, and there are delocalized states for 
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Figure 7. (Color online) The energy dependent density of states is plotted for 
Gaussian distributed random mass for a 30x30 HCL after 10 4 averages for g = 1, m = 2 
(cyan), 1 (blue), 0.5 (red), 0.3 (black), 0.2 (magenta) and (green) in the left panel. 
The right panel visualizes the inverse of the density of states, being proportional to To 
in the variable range hopping model as a function of the energy squared (proportional 
to the carrier density). 

|m| < m c (g) at the NP. The localized states in the tails are described, for instance, by 
the Lippmann-Schwinger equation (1301) . The SCBA with uniform self-energy is not able 
to produce the localized tails. An extension of the SCBA with non-uniform self-energies 
provides localized tails though, as an approximation for large fh has shown [27] . This is 
also in good agreement with our exact diagonalization of finite systems up to 200 x 200 
size. 

A possible interpretation of these results is that there are two different types of 
spectra. In a special realization of m r the tails of the DOS represent localized states. 
On the other hand, the DOS at the NP E = 0, obtained from the SCBA with uniform 
self-energy, comes from extended states [22]. The localized and the delocalized spectrum 
separate at the critical value m c , undergoing an Anderson transition. 

Conductivity: Transport, i.e. the metallic regime, is related to the DOS trough the 
Einstein relation o oc pD, where D is the diffusion coefficient. The latter was found in 
Ref. [22] for E ~ as 

agy/m 2 c -rh 2 2 _ 2 
D = — ^- 0(m c - m ) , (40) 

where a = 1 (a = 2) for MLG (BLG). Together with the DOS po = v/'^9 7r an d the 
scattering rate r] in Eq. (|28|) . the Einstein relation gives us at the NP 

^~0)«A^«£(l-ge (m > -*>)£. (41) 

In the localized regime (i.e. for \m\ > m c ) the conductivity is nonzero only for 
positive temperatures T > 0. Then we can apply the formula for variable-range hopping 
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in Eq. (j2J), which fits well the experimental result in graphane of Ref. [7J- The parameter 
T is related to the DOS at the Fermi level as [Tj 



kBT " K &m • (42) 

where £ is the localization length. T has its maximum at the NP Ep = 0, as shown in 
Fig. [TJ and decreases monotonically with increasing carrier density, as in the experiment 
on graphane [7J. 

In conclusion, we have studied the density of states in MLG and BLG at low energies 
in the presence of a random symmetry-breaking potential. While a uniform symmetry- 
breaking potential opens a uniform gap, a random symmetry-breaking potential also 
creates tails in the density of states. The latter can close the gap again, preventing the 
system to become an insulator at the nodes. However, for a sufficiently large gap the 
tails contain localized states with nonzero density of states. These localized states allow 
the system to conduct at nonzero temperature via variable-range hopping. This result 
is in agreement with recent experimental observations [7J. 
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